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ABSTRACT 

We analyze the gravitational collapse of solids subject to gas drag in a protoplanetary disk. We also 
study the stirring of solids by turbulent fluctuations to determine the velocity dispersion and thickness 
of the midplane particle layer. The usual thresholds for determining gravitational instability in disks, 
Toomre's criterion and/or the Roche density, do not apply. Dissipation of angular momentum allows 
instability at longer wavelengths, lower densities, and higher velocity dispersions than without drag. 
Small solids will slowly leak into axisymmetric rings since initial collapse occurs over many orbits. 
Growth is fastest when particle stopping times are comparable to orbital times. Our analysis of 
particle stirring by turbulence is consistent with previous results for tightly coupled particles, but is 
generalized to loose coupling where epicyclic motions contribute to random velocities. A companion 
paper applies these results to turbulent protoplanetary disks. 

Subject headings: hydrodynamics — instabilities — planetary systems: protoplanetary disks — planets 
and satellites: formation 
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1. INTRODUCTION 

The formation of planetesimals, the first generation of 
solids bigger than a kilometer, remains a poorly under- 
sood, and much debated, stage of planet formation. The 
paradigm of collisional coagulation is very successful in 
explaining the subsequent growth of planetesimals into 
planets. In this stage surfa ce gravity ensures coll isional 
growth of the largest bodies ( Goldreic h et al.l2 QQ4l) . This 
bottom-up "planetesimal hypothesis" tLis sauerl Il993fi 
naturally explains the remnant small bodies and impact 
craters so prevalent in our solar system. Since small dust 
grains coagulate due to electrostatic interactions, it is 
tempting to interpolate that binary collisions are univer- 
sally responsible for the growth of solids. 

However over a wide range of sizes, extending roughly 
from millimeter to kilometer radii, a robust sticking 
mechanism has not been identified, despi te considerable 
exper imental and theor etical effort (s ee lYoudin fc Shi 
120021 hereafter YS02; lYoudinl 12004 and references 
therein). Crossing this gulf of ~ 10 18 orders of magnitude 
in mass is particularly difficult because aerodynamic drag 
induces inward radial drift of sol ids with peak speeds o f 
50 m/s for 0.1 — 1 meter solids ( Wcidcnschil linel 11977ft . 
This would result not only in destructive impacts, but 
also in the loss of solids to the star in ~ 100 yr. 

Collective instabilities involving self-gravity and gas 
drag could directly assemble small solids (perhaps 
millimeter-sized to r esemble known meteoritic inclusions) 
into planetesimals l)Safronovl Il969t iGoldreich fc Ward! 
1973), thereby leapfrogging the regime where collisional 
binding energies are weak and collisional and drift veloci- 
ties are fast. The lone obstacle to gravitational instability 
of solids (hereafter GIS, or GI when referring to gravi- 
tational instability in general) is that turbulence stirs 
solids, lowe ring space densities and increa sing velocity 
dispersions l|WeidenschiHing fc CuzzilH993|) . Turbulence 
is thought to be prevalent i n protoplanetary d isks, mainly 
to drive stellar accretion ijStone et alJl2000|) . Even in a 
passive disk, particle settling generates vertical shear in 



rotational velocities. This can trigger turbulence and im- 
p ede further pa rticle settling ijWeidenschillindl 198Cfl . 

Sekiva (1998) and YS02 showed that turbulence driven 
by vertical shear has a limited ability to loft solids. In 
metal rich disks, with solid to gas ratios several times 
cosmic abundances, not all solids are stirred, resulting 
in midplane GIS. The radial drift mentioned above can 
also serve as a particle enrichment mechanism. Particles 
smaller than 0.1 meter mi grate more slowly and "pi le-up" 
in the inner disk (YS02: lYoudin fc Chianel 12004 here- 
after YC04). These works showed that rapid, dynamical 
GIS is possible. The current work demonstrates that 
slower instabilities exist, and are easier to trigger. 

A common misconception is that GIS only occurs when 
the particle densities exceed the Roche limit or when 
Toomre's parameter, Qt (equation^, is less than unity. 
For single component disks, whether fl uid or particu- 
late, these criteria are well-est ablished l|Toomrel Il964t 
IGoldreich fc Lvnden-Belilll96 5a'). and for thin disks the 
two are usually equivalent (see H2.3(l . Similar criteria 
apply to a dis k of stars and gas, where the coupling is 
gravitational (Rafikov 2001). 

A disk of solids and gas is coupled by dissipative drag 
forces as well. If one assumes perfect coupling between 
solids and gas, no dissipation occurs and the standard 
Roche crit erion for GI is recovered, within facto rs of or- 
der unity l|Sekivalll988l lYamoto fc Sekival I2004D . How- 
ever when slippage is allowed, instability always occurs, 
even if the density is well below the Roche limit. Dissipa- 
tion of angular momentum ensures growth of sufficiently 
long wavelengths . This result will be rederived, but is not 
new ijWardll 19761 120001 . Growth rates are small for low 
densities, but can be much faster than disk lifetimes. The 
mistaken belief in threshold criteria led to the view that 
exceedingly weak turbulence could prevent GIS. Since 
GIS develops at lower densities, stronger turbulence can 
be tolerated. 

This paper investigates the stability of a self- 
gravitating disk of solids (represented as a continuous 
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TABLE 1 
Symbols 



Symbol Reference Eq. 



Meaning 



n 


eq. 01 


orbital frequency 




particle stopping time 


E op ' 




particle surface density 


c, h 


cqs. fWM 


particle random speed, scale height 


Qt, Qr 


cqs. HEJ 


stability parameters for solids 


P 


E//i 


particle space density 


Cg, hg 


eq. <%n 


gas sound speed, scale height 


a 9 


turbulent diffusion parameter 


to,vo,lo 


ro = nto 


eddy turnover time, speed, size 


7 


eq. (TT) 


dimensionless growth rate 


k (fcf) 


wavenumber (fastest growing) 


K (AC f ) 


eq. G3 


dimensionless wavenumber 


Ac 


27r 2 GE/n 2 


traditional GI wavelength 




cylindrical disk radius 






local radial, azimuthal coords. 


Qsonic 


eq. 03 


"mixed" Toomre parameter 



fluid) subject to gas drag in jJH Section estimates the 
particle velocity dispersion and scale height generated 
by turbulent fluctuations. Concluding remarks in ^in- 
clude a summary of previous w orks on GIS in M4.ll A 
companion paper l)Youdinl|2005l hereafter Paper II) ap- 
plies the results of this work to turbulent protoplanetary 
disk models. 

2. DISSIPATIVE GRAVITATIONAL INSTABILITIES 

We investigate the gravitational stability of solids in 
the midplane of a protoplanetary disk as a fluid subject 
to drag aga inst a steady gas backgrou nd. A two fluid 
description ijYoudin fe Goodmanll2005t hereafter YG05) 
would be more general, but our analysis isolates the 
collapse of solids through gas component that is grav- 
itationally stable. Instabilities involving the collapse of 
solids together with gas would be more rapid, but require 
higher densities and are not described here. 

The stability properties of our model are described by 
three dimensionless parameters: t s , Qt, and Qr. Drag 
coupling is measured by 

T s = 17£ stop , (1) 

the particle stopping time, t stopi relative to the orbital 
time 1/17. We will refer to the x s < 1, t s > 1, and 
t s ~ 1 regimes as tight, loose, and marginal coupling, 
respectively. 

Self-gravity is measured by two parameters: 

Qrr- ^ 



Qi 



ttGE 

hn 2 



n 2 



ttGS nGp 



(2) 
(3) 



where £ and p are, respectively, the surface and mid- 
plane space densities of solids, c gives the particle ve- 
locity dispersion, an d h ~ S / p is t he thickness of the 
particle subayer. The iToomrd l)1964f) parameter Qt pits 
the stabilizing effects of pressure and angular momentum 
against self-gravity, while Qr is the ratio of the Roche 
density to the particle density. When c w 17/i, which is 
often the case in astrophysical disks, Qt ~ Qr- In |E1 WC 
show that a distinction is necessary when particles are 
stirred by turbulence. For generality, this section keeps 
t s , Qt, and Qr independent. 



2.1. Basic Equations 

We model the evolution of solids in a localized patch of 
a protoplanetary disk. At cylindrical radius R = R we 
erect a Cartesian coordinate system with x and y axes 
directed in the radial and azimuthal directions, respec- 
tively. The coordinates rotate uniformly at 17 = S7(i? D ), 
the orbital frequency of the solids at R . Our two dimen- 
sional model evolves the surface density and velocity of 
solids: 

£ = E [l + a(x,y,t)] (4) 
V = V (x)+v(x,y,t). (5) 

which are decomposed into steady background, E G and 
V (x), and fluctuating components, where a is a dimcn- 
sioinless overdensity and v = ux + vy. Steady quan- 
tities vary on scales ~ R a (by assumption of smooth 
power-law profiles), and are spatially constant in the 
local approximation. The exception is orbital velocity, 
given to first order in x/R by V D • y — —qtt x, with 
q = —d\n£l/d\nR — 3/2 for Keplerian shear. 

Gas drag induces steady state radial drift of particles 
(inward for radially decreasing gas pressure). In a global 
model, E G would be redistributed in a drift time (YS02, 
YC04), so our local approximation is strictly valid for 
shorter timescales, a condition that will be tested in Pa- 
per II. Since radial motion is spatially constant in the 
local approximation, we remove it via a Galilean trans- 
formation and set V c • x = 0. This is not possible in a 
two-fluid model with different drift speeds (YG05). 

Perturbations evolve via height-integrated continuity 
and Euler equations with a drag acceleration. To lowest 
order in x/ R Q these read: 



dcr dcr 
— - qll x — 
at ay 



V-[(1+«7)V]=0, 



<9v dv „ 

— - qll x- qil uy 

at ay 



217 x v = — V% - 



t'stop 



(G) 
,(7) 



Equations JBJ and Q are analogous t o equa tions (34) 
and (37) in TGoldreich k. Lvnden-Bel l (1965b) where a 
detailed derivation can be found. The linear drag ac- 
celera tion, — v/t s t,nn, enco mpasses Epstein's and Stokes' 
laws ijAdachi et alJ 11976^ and applies to particles with 
radii a < 10(i?/AU) 5 / 4 cm. We dropped small non-linear 
terms C(v 2 ). Our linear and axisymmetric analysis will 
also ignore the C(crv) term and set d/dy = 0. The use of 
fluid equations is valid for processes involving timescales 
longer than t stop - Thus loosely coupled solids can some- 
times be t reated as a fluid, for instance in calculating 
drift rates l)Nakagawa et alJ[l986|) . 

The effective potential, x = & + n, includes the per- 
turbed gravitational potential, and an effective pres- 
sure (more precisely, enthalpy): 



n = c 2 ln(l + a) 



(8) 



The approximation that the velocity dispersion of solids 
acts as a thermodynamic pressure is a useful, if inex- 
act, analogy that has bee n explored extensively in colli- 
sionless stellar dynamics ijBinnev fc Tremainelll987T) . Its 
appropriateness for the current problem should be in- 
vestigated in more detail. Our choice of an isothermal 
equation of state is not crucial. A different choice, e.g. 
II = c 2 (l + cr) n , would give similar results for linear per- 
turbations. 
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The gravitational potential of an individual Fourier 
mode is: 



$ = -2-KGT, -T(kh) 



(9) 



where f(k x , k y ) is the 2D Fourier transform of f(x, y) and 

k = ^Jkl + k$. The softening term, T(kh) = 1/(1 + kh) 

mimics finite disk thickn e ss in a ve rtically integrated 
model (|Vandervoortlll97ft IShul ll984). This factor con- 
nects the thin disk limit for long waves, fc/i < 1, to the 
three dimensional result for fc/i> 1. 

2.2. Axisymmetric Dispersion Relation 

We find linear, axisymmetric solutions to equations 10 
El by assigning a Fourier dependence, exp[i(kx — ujt)], 
to perturbed quantities. The resulting dispersion rela- 
tion reads: 



(uj + i/t Btop ) [V 2 - 2irGEkT(kh) + k 2 c 2 
-uj (uj + i/t stop )] - itf/tstop = , 



(10) 



where above and henceforth we drop the naught sub- 
scripts on background quantities. Equ ation fliup ag rees 
with the di spersion rela tions derived in lWardl l)1976l eq. 
A-22 ) and I War dl l|2000L eq. 31). We introduce a dimen- 
sionless growth rate and wavenumber, 

(11) 

nGZk/n 2 , (12) 

where is the characteristic gravitational wavelength, 
and express the dispersion relation as a real cubic: 



7= —luj/Q 
IC = kX G /(2n) 



where 



F(K,Qt,Qr) 



1 



2/C 



1 + KQ R 



T7 1 = 0, 



QtK. 2 



(13) 



(14) 



The function F encapsulates all relevant physics (rota- 
tion, self-gravity, and velocity dispersion) except for dis- 
sipation. Since dj/dF = (dF/dj)^ 1 < for 7 > 0, 
growth is faster for lower values of F, i.e. stronger self- 
gravity. Note that the minimum of F(JC) is always < 1, 
a fact we will exploit shortly. 

2.3. Non-dissipative GI 

We begin by showing that our model reproduces stan- 
dard non-dissipative results. As t s — > 00, equation IJ13JI 
reduces to 7 2 = — F plus a neutral mode with 7 = uj = 0. 
Thus instability requires F < 0. A more familiar, but 
identical, form for the density wave dispersion relation 
reads: 



= fr — 2irGmT(kh) + k 2 



(15) 



In the thin disk limit T — > 1, or equivalently Qr — » 0, we 
recover the familiar instability criterion Qt < Qt, cmt — 
1, and the marginally stable wavenumber is /C cl -it = 1. 

It is worth considering finite thickness effects in the 
dissipation-free case. In a cold disk with Qt = 0, sta- 
bility is assured when Qr > 2. The Roche limit is en- 
forced in the absence of dissipation. For less extreme 
cases, as Qr increases (from — ► 2), <5r,crit decreases 
(1 — > 0), and /C cr ; t increases (1 — > 00). Physically, at 
lower volume densities the disk must be colder to go un- 
stable, and it will do so at shorter wavelengths. If we 



require Qji = Qt (the standard c = £lh relation) then 
Qr.crit ~ 0.55 and /C cl -it ~ 1.2. This actually overesti- 
mates the stabilizing influence of finite thickness. Self- 
gravity shrinks the thickness of an isothermal fluid so 
that Q T = Q R y/l + 2/Q R i|Wardll200ft our eq. [57j) in- 
cludes the same correction for vertical self-gravity). This 
gives Qr.crit ~ 0.77 and /C cr it ~ 1.05, a less dramatic 
finite thickness correction. 




10-31Q-2 0.1 1 10 10 2 10 3 



i cncy ) 

for dissipative GI as a function of stopping time and F (eq. 1141 1. 
The tight-c ouplin g approximation is overplotted with dotted con- 
tours. See H2.4.1I 



2.4. Dissipative GI 
2.4.1. Dissipation Guarantees Instability 

Drag forces damp angular momentum that would oth- 
erwise support long wavelength modes. This qualita- 
tively changes stability characteristics. While F < 
was needed for instability in the absence of dissipation, 
equation H13JI with r s real and positive implies growth for 
F < 1. Appendix 1X1 proves this necessary and sufficient 
condition. 1 Since F < 1 as JC — > (eq. |14|). sufficiently 
long wavelengths are always unstable. 

Figure ^ shows growth rate contours as a function of 
F and t s . The stability boundary at F = 1 is clear. At 
large r s , contours of slow but finite growth (7 < 0.1) 
converge to F = 0, the dissipation-free result. Growth 
rates peak around t s — 1 because dissipation is strong 
enough to damp vorticity on an orbital time, but not so 
strong that it impedes collapse. 

2.4.2. Fastest Growing Wavenumbers 

The wavenumber giving the fastest linear growth rate, 
k{ (or the dimensionless /Cf), is the most relevant until 

1 Appendix ^ also shows that only one of the three roots can 
grow, and that this mode does not oscillate. 



4 



A.N. Youdin 



non- linear effects become significant. Since all K, depen- 
dence arises through F, JC{ is the root of: 



d'f d'f dF 

dIC = dFdJC = ' 



(16) 



Since dj/dF < for all growing modes, K,f is just the 
root of dF/dJC = which satisfies: 



K,fQ T (l+K t Q R ) 2 = 1. 



(17) 



Thus /Cf is determined by the stability parameters, seem- 
ingly independent of dissipation. Without dissipation the 
longer waves would not grow. Also drag coupling influ- 
ences Qt and Qr values (jJSJ. 

Modes will be long (or short) compared to the disk 
scale height if kfh — ICfQa < 1 (or > 1). We refer to 
these as the thin (or thick) disk limits, respectively. Note 
that kfh depends on = Q t /Qr only. Series solutions 
show: 

§-^ +0 ( Q ~ 3 ) ife>>1 ' ( 18a ) 



kfh 



ev 



i--l + o(e^ 3 ) ife«i. (i8b) 



The -1 / 3 dependence in equation (|18b(l means that ex- 
treme conditions (very cold and low density) are needed 
for waves to be much shorter than h. 

The fastest growing wavelengths, Af = 2n/kf = Ac//Cf, 
satisfy: 



Af_ 



(QtQr) 



2/3 



if kfh < 1, 
if kfh > 1. 



(19a) 
(19b) 



The thin disk limit is independent of scale height as it 
must be. In both regimes, hotter disks (Qt large) disks 
go unstable at longer wavelengths to avoid "pressure" 
stabilization. Most importantly, modes will be much 
larger than Xq when the stability parameters are large. 

2.4.3. Tight Coupling 

Particle sizes and gas densities relevant to planetesimal 
formation usually correspond to tight coupling, t s -C 1. 
In this regime a destabilized neutral mode grows at a 

7 ~r.(l-F)+0(7f). (20) 
The two other modes are less interesting: strongly 



damped inertial oscillations with 7 



-t7 1 ±i- The dot- 



ted contours in Figure ^ show that equation l|20(l agrees 
with the full solutions for small t s . Leading order solu- 
tions for ICf (equations m equations 1)141120(1 give the 
fastest growth rates as: 



7r 



—p if kfh <C 1 , 
Qt 



2ts_ 
Qr 



if kfh > 1. 



(21a) 
(21b) 



Growth in the thin disk limit depends on the velocity 
dispersion but not the scale height, and vice- versa in the 
thick disk limit. For large values of the stability pa- 
rameters, dissipative growth is slower than the particle 
settling rate, CIt s . Thus linear growth proceeds through 
a series of quasi-equilbrium states. The loose coupling 
case is evaluated in Appendix iBl 



2.5. Why Dissipation Helps 

To explain how angular momentum dissipation aids in- 
stability, we offer the following analogy with the Rayleigh 
instability. The arguments are fairly technical and can 
be skipped if desired. Generating axisymmetric density 
enhancements requires radial velocity perturbations, u. 
In the absence of dissipation, the azimuthal force equa- 
tion conserves angular momentum: 



dv/dt= (-2 + q)ttu 



(22) 



or v = — (2 — q)u/j, in terms of the dimensionlcss growth 
rate. The Coriolis force is partly blunted by differential 
rotation, with q — 3/2 in a Kepler disk. For slow growth 
(7 <C 1), angular momentum conservation requires \v\ 3> 
|u| unless q « 2. Azimuthal motions in turn supply a 
radial Coriolis force, fc, r — 2^w = —2(2 — q)uil/j. For 
q < 2, fc, r decelerates u and has the form of a drag force 
that becomes stronger for slower growth. For q > 2, 
the sign of the radial acceleration changes, giving rise 
to the well-known Rayleigh instability for an outwardly 
decreasing angular momentum gradient. 

Now we include actual drag forces and restrict our- 
selves to q = 3/2. Azimuthal force balance now includes 
a loss term for orbital angular momentum. 



dv/dt = — flu/ 2 — v/tgtop. 



(23) 



If t s 7 <C 1, i.e. strong drag or slow growth, the Corio- 
lis force (blunted by differential rotation) is balanced by 
drag. The resulting azimuthal perturbation, v « —t s u/2, 
is much smaller in magnitude than the drag free case 
(where v ~ — u/7) as if q were approaching 2. In turn 
fcr = ~ —£It s ii (the prime differentiates between 
this and the drag-free case) offers much less resistance 
to collapse. For r s <C 1, the radial component of the 
drag force, fd, r — —u/ts, is the dominant obstacle to 
collapse. Since \fd, r \ <C |/c>| ~ M^/7 f° r 7 ^ T s, the 
presence of drag allows slow growth that would otherwise 
be prevented by angular momentum conservation. Since 
this collapse is primarily resisted by drag, it strongly re- 
sembles particle settling (to the center of the contracting 
annulus) at the terminal velocity. 

3. PARTICLE STIRRING BY TURBULENCE 

The previous section treats Qt and Qr as free parame- 
ters. To constrain these stability parameters, this section 
calculates the velocity dispersion, c, and scale height, h, 
of solids coupled to turbulent fluctuations. Our approx- 
imate expressions reproduce known results in the tight 
coupling limit, but we generalize the analysis to arbitrary 
stopping times. In contrast to most astrophysical disks, 
we find c 7^ ilh, with c <C £lh for weak turbulence and 
tight coupling. This motivated our use of two stability 
parameters. We ignore other contributions to random 
motion, like gravitational scattering, aerodynamic drift, 
or phyical collisions, which are less significant, particu- 
larly for small, tightly coupled solids, see jjn] 

3.1. Turbulent Diffusion 

We assume that gaseous turbulence acts a diffusive vis- 
cosity 



/n 



0^0 ■ 



(24) 



defined in terms of a 9 , a dimensionlcss transport coef- 
ficient for gas, and c g , the sound speed. The second 
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equality relates v g to the characteristic eddy speed, vq, 
and turnover time, to- We must specify two of a g , vq, 
and/or to to describe the turbulence. 

In general, particl es are not per f ectly coupled to tur- 
bulent fluctuations. iMizuno et alJ 1(19881) calculated the 
RMS velocity of solids, c„, in response to a Kolmogorov 
spectrum of eddies. They found c„ ~ vq for parti- 
cles well coupled to the largest eddies (i s top < to), and 
c v ps wo \/to/tstop if t s top > to ■ Both regimes are covered 
by: 



t'o 



y/l + T s /t 



TO + T s 



(25) 



where To = flto- 

The diffusivity of the solids, v s = a s c 2 /Vl 1 is defined 
analogously to v g . Assuming the velocity perturbations 
(equation I25|) fluctuate on the forcing time, to, we have 
v s w c^tQ. The Schmid t number, a ratio of particle to 
gas diffusivities, is then ijCuzzi 



Sc = v g /v s = a g /a s « 1 + t stop /t 



o ■ 



(26) 



As i st0 p - ¥ 0, solids diffuse as efficiently as gas in this 
approximation. 

3.2. Particle Scale Height 

We estimate an equilibrium scale height for the solids 
by equating diffusion times, tdiff = h 2 /v s with settling 
times, tsett< Tightly coupled particles fall to the midplanc 
at terminal velocity and f se tt = t s top) j where self- 



gravity amplifies vertical gravity by a factor 
V> = 1 + 2/Q R . 



(27) 



For weak GI with Qr 3> 1, we will frequently set ip = 
1. Loosely coupled particles undergo damped vertical 
oscillations with t se tt = istop, with no correction for self- 
gravity. The combined result for arbitrary f st0 p is: 



tsett — t. 



stop 



l/(^ 2 tsto P ) • (28) 



Equating tdiff = tsett and using equation 1(26(1 gives h, 
here expressed relative to the gas scale height, h g = c g /f2: 



h 
h 



1 + ^{h) T 2 



'i>(h)T s (l+T s /T ) 



(29) 



Our diffusion time implicitly assumed that h < h g . If 
equation l|29|l gives h > h g , then particles are evenly 



mixed with gas and h 



The h dependence of ip 



is noted, because when self-gravity is included, equation 
(|29|) defines h implicitly as the root of a cubic. In all 
cases, h increases with a g . 

As t s — > co, h — > h g y/a g To — Voto, a constant equal 
to the characteristic eddy size, Iq. One might naively ex- 
pect h to decrease at large t s , but in our analysis both 
the settling time and the diffusion time increase cx t s . A 
non-linear contribution to Sc would drive h — * in the 
loose coupling limit. Since h ~ Iq for r s 3> 1, the approxi- 
mation that turbulence acts diffusively is only marginally 
justified in the loose coupling regime. For tight coupling, 
h 3> Iq, and the diffusion approximation is well justified. 



3.3. Velocity Dispersion 

The random motions of the solids includes direct kicks 
from turbulent fluctuations, c„ (eauation l25|) . plus a con- 
tribution from epicyclic oscillations, c orD . We assume 
that eccentricities, e, and inclinations, i, are small and 
in the dispersion dominated regime, e ~ i <C 1. The net 
velocity dispersion, 



u orb 



(30) 



is dominated by c v for tight coupling, and c or b for loose 
coupling. 

To estimate c or bj consider epicyclic oscillations of an 
individual particle with gas drag : 

x = — Q 2 x — x/ts , (31) 

We isolate radial motions since we are concerned with 
support against radial collapse fsee iGaraud et all 120041 
for a detailed treatment of vertical oscillations). The 
epicyclic frequency needs no correction for self-gravity 
since the total disk mass is small compared to the central 
star. Equation 1)31(1 describes a damped harmonic oscil- 
lator. Driving by turbulent fluctuations is not included, 
but should maintain a constant amplitude h w iR w eR. 2 
We define c or b = w max , where v m&x is the maximum 
speed achieved by a particle released at rest from x — h. 
For t s ~> 1, "Umax ~ f2/i as the particle crosses x = 0, 
the radius of the orbit's guiding center. With t s <C 1 the 
particle reaches terminal speed w ma x ~ Qr s h shortly after 
release. Intermediate cases reach a maximum speed: 



Corb 



flhr s 

1 + Tt 



4> 



TO + T s 



1 



(32) 



using equation l|29|l to eliminate h. 

Comparison of equations 125fl and (|32|1 shows that c or b 
dominates random motions for r s > 1/ro, while c„ is 
larger for tighter coupling. The dynamically interesting 
ratio 



c 

Qh 



Qr 



^{l+T s f 



(33) 



I -!- r : , V ' t t s {\+Vt 2 ) 

If particles are well coupled to eddies, r s < to, and to 
orbital motions, r s <C 1, then c/(0/i) w \/iPt s /t -C 
1. In the loose coupling limit we recover the standard 
c/(Clh) w 1. For eddies with short turnover times, To <C 
1, c/(flh) w y2/ro > 1 at t s = 1, showing that a wide 
variety of behaviors are possible. 

3.4. Results: Slow vs. Fast Eddies 

Turbulence cannot be described by a single free pa- 
rameter, a g . We must also specify the speed of eddies. 
We consider two cases, though intermediate speeds are 
possible. 

3.4.1. Orbital Turnover Times 

If orbital shear is involved in the generation or destruc- 
tion of turbulence one naturally expects To = 1. Since 
eddies with longer turnover times are destroyed by ra- 
dial shear, this gives the slowest allowed eddy velocity, 

2 Ideally turbulent driving and orbital dynamics would be 
treated concurrently, a task beyond the scope of the current work. 
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Fig. 2. — Particle scale height, h (solid curves), and veloc- 
ity dispersion, c (dashed curves), vs. stopping time for turbulence 
with orbital turnover times (black curves) and faster eddies (grey 
curves) with ct g = 10~ 4 in both cases. For orbital turnover times, 
the two contributions to c are shown separately. With faster eddies, 
c v dominates and only the total c is plotted for clarity. 



v a = Jc7 g ~c g . Weak turbulence with a s <C 1 is very sub- 



sonic. Values for h and c, plotted in Figure [21 with black 
curves, are obtained from equations 1|29[) . (|25l) . and l|32[l 
with To = 1 and if> = 1. When a g < r s < 1, coupling is 
imperfect enough to allow some settling and we recover 
the standard result (Dubrullc ct al. 1995, also in YC04): 



h 



(34) 



For loose coupling h w ^u g ~h g . The velocity dispersion 
maintains a nearly constant value c ~ ^JcT g ~c g , but the 
source shifts from direct kicks by turbulent fluctuations 
for t s <C 1 to epicyclic motion for t s 1. For tight 
coupling c/(Q,h) w -^77 <C 1 because small velocity kicks 
loft particles efficiently. 

Expressions for the stability parameters are simplified 
by defining 

,. ( ) 

(35) 



Cgfl 

ttGS' 



a "mixed" Toomre parameter involving the gas sound 
speed and particle surface density.The stability parame- 
ters are approximated by 



if t s 
if a. 



< a g 

< T s < 1 



if T S > 1 



(36a) 



(36b) 



Numerical solutions are plotted for several values of a g 
in Figure [3] (top). Stability parameters are larger, i.e. 
self-gravity is weaker, for stronger turbulence. For tight 
coupling we confirm Qt/Qr — c/(Slh) ~ <C 1. 
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Fig. 3. — Top: Stability parameters Qr (black curves) and 
Qt (grey curves) vs. stopping time for orbital turnover times with 
Qsonic = 5 X 10 3 . The gas diffusivity, ct g , takes values of 10 — 8 
(solid curves), 10 — 6 (dashed), and 10 — 4 (dot-dashed). The dotted 
line orients the reference value of 1. Bottom: same as above but 
for faster eddies. 



3.4.2. Fast Eddies 

We now consider what happens if turnover times are 
shorter and eddy speeds are faster. We still want eddies 



to be subsonic for a g <C 1, so we choose i>o 
this case turnover times are short, tq = ,/i 



I/ 4 



In 



<C 1, and 



eddies are small, Iq ~ a g h g . For comparison recall that 
with orbital turnover times we had tq — 1, vq = Jouc 



and Iq — ^Ja g ~c g . We favor the case of orbital turnover 
times, and the case of fast eddies (FE) is primarily meant 
to be instructive. 3 

Figure |21 (grey lines) plots c and h for FE. Compared 
to turbulence with orbital turnover times (black lines) 
h is smaller, and thus densities higher, because particles 



are less tightly coupled to eddies once r s > r 

3 We also investigated the extreme case of vq = c g , and the 
results were qualitatively similar to FE. Since sonic turbulence im- 
plies shocks, our formalism might not be appropriate, and we do 
not present those results. 
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For t s < 1 FE generate a larger velocity dispersion. We 
know from i j2.4.3l that, in the long wave limit, larger c 
gives slower growth even if h is smaller. Paper II confirms 
that this is the correct regime to consider. Thus growth 
is slower for tightly coupled particles for FE. 4 For t s > 1 
however, particles are so poorly coupled to FE that c is 
lower despite the higher eddy speed. Thus for loosely 
coupled solids GIS is more rapid for FE (with a g fixed). 

Figure|31( bottom) plots the stability parameters for the 
fast eddy case. For weak turbulence and loose coupling, 

values of Qr m ag^ 4 Q S onic are remarkably small. How- 
ever, as long as Qt ^ 1 violent dynamical instability is 
prevented. Both plots (Fig. [SJand Fig. |3J) show c 3> (fi/i) 
for fast eddies in the range ^fa^ < t s < 1/ y/ot^. Thus 
FE produce a thin, hot disk. 

4. DISCUSSION 

4.1. Previous Work 

Ma ny authors hay e con sidered GI of solids in a gas 
disk. iSafronovl l)1969j) and iGoldreich fe Wardl l)1973fl are 
most associated with planetesimal formation by GI. Nei- 
ther included dissipation by drag in their stability crite- 
ria. Instead they used the Roche and Toomre criteria, 
respectively. Both of these give GI at the characteris- 
tic wavelength Aq. In minimum mass models this leads 
to the kilometer estimate for planetesimal sizes. The 
use of threshold criteria led later researchers to conclude 
that turbulen ce could readily preve n t planetesimal for- 
mation by GI HWeidenschillinell 98(1). IGoldreich fc WaTdl 
(1973) did use drag to explain dissipation of angular mo- 
mentum in the second stage of their two-stage collapse 
model. Moreover, they noted that "gas drag does desta- 
bilize axisymmetric perturbations for wavelengths larger 
than (Ag)," but did not investigate the possibility in de- 
tail. 

Subsequent work has included dissipation to find GI at 
low densities, but the impli cations have not been widely 
appreciated. iSpieeell (^972) considered the Jeans insta- 
bility, i.e. gravitational collapse in a uniform non-rotating 
medium, in a coupled two fluid system. Spiegel de- 
scribed " clumping of dust a lone" as it "slips through 
the gas." IWardl l|197fit 1200(1) considered the effect of a 
linear friction term on GI in disks. He found "destabi- 
lized neutral modes" always exist at long wavelengths. 
Our fJ2 is a more detailed investigation of the same ba- 
sic model. Ward analyzed the dispersion relation in the 
long- wavelength limit, which allowed him to avoid esti- 
mating the velocity dispersion or scale height. Since we 
isolate the fastest growing mode, our growth rates are 
s omewhat faster. 

iCoradini et alJ l)1981|) investigated GI in a disk of two 
coupled fluids, but the results are confusing. A minimum 
density threshold for GI is claimed, even though they 
derive a growth time that reproduces the IWardl l)197fif) 
result with no density cutoff. The problem appears to 

4 This conclusion may not hold when collisional dissipation, 
which Appendix |C] shows is significant for c 2> Oft, is included. 



be that ICoradini et all (|198 If) arbitrarily impose a maxi- 
mum wavelength cutoff by appealing to the "conservation 
of angular momentum of the grains" even while correctly 
noting that "angular momentum can be freely exchanged 
w ith the stationary ga s." 

iGoodman fc Pindorl l)2000[) primarily considered drag 
instabilities. They also found long wavelength GI with 
slow growth rates for certain choices of nebular parame- 
ters. Since they use a different drag prescription (bound- 
ary layer or "plate" drag) and include diffusive effects, 
direct comparison with our work is difficult. 

The response of solids to a gravitationally u nstable 
gas disk has been considered. iNoh et all l|199lD found 
that adding solids inc reased the growth rate of global, 
non-axisymmetric GI. iRice et all (j2004f) showed that 1- 
10 meter planetesimals accumulate in the spiral arms of 
a massive gas disk. The literature also contains many 
mechanisms for concentrating solids without self-gravity 
(see YS02 for an incomplete summary). 

4.2. Summary and Future Work 

This work investigated the gravitational collapse of 
solids subject to gas drag. We explored the dependence 
on layer thickness and velocity dispersion (normalized 
to the strength of self-gravity by Qr and Qt, respec- 
tively), and the aerodynamic stopping time, r s . Even 
when self-gravity is weak, sufficiently long wavelengths 
are always unstable to collapse over many orbital times. 
This initially forms a large circumstellar annulus which 
should subsequently fragment to form less massive plan- 
etesimals. 

We also investigated particle stirring by turbulence, 
which should largely determine the velocity dispersion, 
c, and scale height, h. In contrast to the usual c ~ Qh 
relation, we found c <C flh for tightly coupled particles 
when eddy turnover times are orbital. Paper II uses these 
results to investigate the gravitational accumulation of 
solids in turbulent disks. 

Our dynamical model could be refined by adding 
physics such as a dispersion of particle sizes, a two-fluid 
treatment where the gas responds dynamically, general- 
ization to three dimensions, and non-linear growth. Stir- 
ring by turbulent fluctuations can also be analyzed in 
more detail, including the use of direct numerical simu- 
lations. 
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APPENDIX 

A. PROOF OF STABILITY CRITERION 



Proof of the F < 1 stability criterion follows. We express the dispersion relation as 

P( 7 ) = 7 3 + 2 7 2 /r s + (t- 2 + F) 7 + (F — l)/r s = . 



(Al) 
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Fig. A4. — Growth rates (top) and oscillation frequencies (bottom) of growing (solid curve) and damped (dotted curves) modes in the 
loose coupling regime. With F = —0.1 < 0, self-gravity is strong enough to produce GI without dissipation. Bifurcation corresponds to 
the changing character of modes. See JE] 



Properties of the real roots follow readily from Descartes' "rule of signs" (J. Goodman, personal communication). 
With F < 1 there is one sign change to the coeffients of P(^f) and thus one and only one purely growing mode. For 
F > 1 there are no purely growing modes, and the coefficients of P(— 7) show that there are 1 or 3 purely damped 
modes. 

To complete the proof we show that there can be no complex roots with a positive real part, which is true for all 
F. Suppose that P has one real root, 71, and a pair of complex conjugate roots, 72 and 73 = 7^. From (7 — 71) (7 — 
7a) (7 ~ 7a) = -f (7)) the quadratic term requires that the real part of the complex roots ^(72) = ^(73) = — 7i/2 — T.T 1 . 
We (fallaciously) suppose that ^(72) > 0. Then 71 < — 2/r s , requiring F > 1 and P(— 2/r s ) > (with only one real 
root, P > for all 7 > 71). Since P(—2/t s ) = — (1 + F)/t s < 0, the contradiction proves that there are no 

overstable oscillations (growing complex roots). Thus F < 1 is a necessary and sufficient condition for instability, and 
growing modes have no oscillatory component. 

B. LOOSE COUPLING 

The loose coupling limit may be relevant to planetesimal formation in the outer disk or for the subsequent assemblage 
of larger bodies like planetesimals or planetesimal fragments. It also establishes the link between destabilized neutral 
modes and unstable density waves. If t gIOW = 1/(7^) < t stop , i.e. 7T S > 1, the fluid approximation is poor and a 
kinetic theory treatment would be more appropriate. When t s 3> 1 approximate formulae for the growth rate fall in 
one of three regimes: 




l + F 
-2Ft x 



Ft, 



if F<-r- 2/3 , 
if r7 2/3 < F, 



-1/3 



FtI/ 3 /S if \F\ < r7 2/3 . 



(Bla) 

(Bib) 
(Blc) 
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The first case, with F < but not too close to 0, is an unstable density wave with a small drag correction that increases 
the growth rate. This is seen in Figure H where contours are nearly horizontal for F < 0, t s ^S> 1 . For case two, F > 
by a finite amount, and the neutral mode is destabilized for F < 1. Case three shows that the transition across F = 
is well-behaved. 

2/3 | _ | 

This transition is associated with a bifurcation point at F ~ — 0.52r s , see Figure lA41 To the left of the bifurcation 
point, as r s — > oo, the unstable mode is an ordinary density wave. To the right of the bifurcation point, the unstable 
mode is characterized as a neutral mode, and for r s < 1 connects with equation (|2U[1 . Note that the bifurcation 
actually involves the damped modes. When F > 0, there are no unstable density waves, or bifurcation points, and the 
unstable mode can always be identified as a neutral mode. 

Growth rates of the fastest growing modes in the F > regime are approximated by 



1 



7f 



T S Q T 

2 



t s {Qi 



if jfef/K 1, (B2a) 
if kth » 1, (B2b) 



where Qr > 2. As in the tight coupling case, growth in the long (or short) wave regime depends primarily on Qt (or 
Qr), respectively. 

C. VELOCITY DISPERSION: IGNORED EFFECTS 
Aerodynamic Drift 

Since we have used a single particle size, aerodynamic drift is uniform and has no effect on random speeds. General- 
ization to a dispersion of particle sizes (or more precisely, stopping times) would give a range of drift speeds, with the 
radial component being relevant for axisymmetric collapse. We now show that drift speeds are smaller than velocities 
induced by turbulent forcing, c, except per haps near marginal cou pling. 

The inward radial drift speed of solids is l)Nakagawa et al.lll98fj YG05): 

V ^ = ^Pg o T 2 ( C1 ) 

where p g = p g /(p + p g ) is the midplane gas fraction and r] ~ Cg/vI c measures radial pressure support. We ignore 
the effect of turbulent diffusion on radial flow ijTakeuchi fc Linll2002|) . For t s -C 1 and c » yfa^c g f H3.4.1fl a bit of 
manipulation shows: 

c _Wcj h v K /c g E v K /c g 

Where we used equation and \i v = p/ptot = 1 — Pg is the solid mass fraction. The inequality holds because 
even with S/S g ~ .01 in the final expression, the contributions to the second fraction are all > 1 and will dominate. 
Surprisingly, turbulent speeds dominate even for very weak turbulence. This is because the particle layer becomes so 
thin and dense that i^r decreases even faster due to the p 2 g term. For r s ^> 1, we can similarly show that radial drift 
is ignorable, but drift could contribute to the velocity dispersion of marginally coupled solids. 

Collisional vs. Aerodynamic Dissipation 

We have ignored the role of inelastic collisions on the velocity dispersion, an effect that could make the di sk more grav- 
itationally unstable. The re lative importance of collisions and drag depends on the applicable drag law ijAdachi et alJ 
119761 IWeidenschillindll977j) . the relative abundances of solids and gas, and the ratio c/(Cth). In the Epstein regime, 
the ratio of drag to collisional timescales, 

!l°R _ PsO_J^c_ _ Ji_ < 1 / C3 s 
tcoil PgCg hp s a Y.gVlh 

is independent of particle size, a, and drag dominates when gas is more abundant than solids (S 9 > S) and/or c -C 
as found in t|3.4.1l If turnover times are short f ij3.4.2|) and c ^> ilh then inelastic collisions would provide additional 
dissipation. 

In the Stokes regime, particles are bigger than the gas mean free path, a > A m f p , but not large enough to trigger 
turbulent wakes. In this case collisions are relatively more important: 



4a 



j.St y 
■-stop t-' 

icon Eg Qh 9A m f p 



(C4) 



However increasing a likely brings particles to the turbulent drag regime before collisions dominate. For large particles 
in the turbulent drag regime: 

t stop _ 8 S c c g ^ ^ S v K ^ C5 ^ 
icon 3Cd Eg flh Av E ff c g 
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In the second equality we use the drag coefficient Cd — 0.44, and since t s 3> 1 for turbulent drag, 5 we take c ss flh and 
the relative velocity Av = t}Vk ~ c^/vk- Thus for turbulent drag, collisions and drag have comparable importance at 
cosmic abundances £ s w 100S, and collisions are more significant when gas is depleted. 

Our dynamical model uses a linear drag law, either Epstein's or Stokes'. In these regimes collisional dissipation is 
significant if turbulent stirring is vigorous and c 3> £lh. Thus solids may be more gravitationally unstable than our 
model indicates. Of course when particle overdensitics arc highly non- linear (a stage we do not describe in this work) 
collisional damping will dominate drag. 
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